home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Environments / MacMETH 3.2.4 / More Examples / Newton.MOD < prev    next >
Text File  |  1996-06-20  |  6KB  |  242 lines
  1. MODULE Newton; (* HS 19.10.91 *)
  2.  
  3.   (* NOTE: This module can only be compiled with the compiler "Compile20" *)
  4.  
  5.   FROM SYSTEM IMPORT Exp, Ln, Sqrt, Sin, Cos, ArcTan;
  6.   FROM Terminal IMPORT BusyRead;
  7.   FROM CursorMouse IMPORT GetMouse;
  8.   FROM InOut IMPORT Read, Write, WriteString, WriteLn;
  9.   FROM Windows IMPORT SetWindow, ResetWindow;
  10.   FROM GraphicWindows IMPORT Window, OpenGraphicWindow, CloseGraphicWindow, Clear,
  11.                              SetPen, MoveTo, Dot, IdentifyPos;
  12.  
  13.  
  14.   CONST NrPixels = 256;
  15.  
  16.         (* colours *)
  17.         white   = 449;
  18.         yellow  = 65;
  19.         green   = 321;
  20.         zyan    = 257;
  21.         blue    = 385;
  22.         magenta = 129;
  23.         red     = 193;
  24.         black   = 1;
  25.  
  26.  
  27.   VAR ch             : CHAR;
  28.       v,w            : Window;
  29.       environ        : BITSET;
  30.       ox,oy,nx,ny    : INTEGER;
  31.       dx,dy,max,p,q  : REAL;
  32.       xmax,xmin      : REAL;
  33.       ymax,ymin      : REAL;
  34.       ux,uy,uux,uuy  : REAL;
  35.       deltax         : INTEGER;
  36.       maxiter        : INTEGER;
  37.       this           : INTEGER;
  38.  
  39.  
  40.   PROCEDURE f(x: REAL): REAL;
  41.   BEGIN
  42.     RETURN x * Cos(6.0*x) * Exp(-x*x)
  43.   END f;
  44.  
  45.   PROCEDURE ForeColour(c : LONGINT); CODE 0A862H;
  46.   PROCEDURE BackColour(c : LONGINT); CODE 0A863H;
  47.  
  48.   PROCEDURE Colour(k : INTEGER);
  49.   VAR c : LONGINT;
  50.   BEGIN
  51.     IF k >= maxiter THEN c := black
  52.     ELSE
  53.       k := k MOD 16;
  54.       CASE k OF
  55.           0 : c := white;
  56.         | 1 : c := yellow;
  57.         | 2 : c := yellow;
  58.         | 3 : c := green;
  59.         | 4 : c := green;
  60.         | 5 : c := zyan;
  61.         | 6 : c := zyan;
  62.         | 7 : c := blue;
  63.         | 8 : c := blue;
  64.         | 9 : c := blue;
  65.         |10 : c := magenta;
  66.         |11 : c := magenta;
  67.         |12 : c := magenta;
  68.         |13 : c := red;
  69.         |14 : c := red;
  70.         |15 : c := red;
  71.       END;
  72.     END;
  73.     ForeColour(c);
  74.   END Colour;
  75.  
  76.   PROCEDURE GetPos(VAR i,j: INTEGER);
  77.   CONST ML = 15;
  78.   VAR mouse : BITSET; x,y : INTEGER;
  79.   BEGIN mouse := {};
  80.     REPEAT GetMouse(mouse,x,y) UNTIL NOT(ML IN mouse);
  81.     REPEAT GetMouse(mouse,x,y) UNTIL ML IN mouse;
  82.     REPEAT GetMouse(mouse,x,y) UNTIL NOT(ML IN mouse);
  83.     IdentifyPos(w,x,y);
  84.     i := x; j := y;
  85.   END GetPos;
  86.  
  87.   PROCEDURE SANE(VAR e: BITSET; OpWord: CARDINAL); CODE 0A9EBH;
  88.  
  89.   PROCEDURE SaveFPEnv;
  90.   BEGIN SANE(environ, 3)
  91.   END SaveFPEnv;
  92.  
  93.   PROCEDURE ClearFPEnv;
  94.   VAR e: BITSET;
  95.   BEGIN e := {}; SANE(e, 1)
  96.   END ClearFPEnv;
  97.  
  98.   PROCEDURE RestoreFPEnv;
  99.   VAR e: BITSET;
  100.   BEGIN e := environ; SANE(e, 1)
  101.   END RestoreFPEnv;
  102.  
  103.  
  104.   PROCEDURE Calculate;
  105.   VAR x,y: REAL; ix,iy: INTEGER;
  106.  
  107.     PROCEDURE iteration(pe,qe,xe,ye: REAL) : INTEGER;
  108.     VAR a,y1,y2,x1,x2: REAL; counter: INTEGER;
  109.  
  110.     BEGIN
  111.       IF xe = ye THEN RETURN maxiter END;
  112.       x1 := xe;
  113.       x2 := ye;
  114.       counter := 0;
  115.       y1 := f(x1);
  116.       y2 := f(x2);
  117.       REPEAT
  118.         a := y1 - y2;
  119.         IF a = 0.0 THEN RETURN maxiter END;
  120.         a := a / (x1 - x2);
  121.         x2 := x1; y2 := y1;
  122.         x1 := x1 - y1/a;
  123.         y1 := f(x1);
  124.         INC(counter);
  125.       UNTIL (ABS(y1) < max) OR (counter >= maxiter);
  126.       RETURN counter
  127.     END iteration;
  128.  
  129.   BEGIN (* Calculate *)
  130.     SaveFPEnv;
  131.     ClearFPEnv;
  132.     dx := (xmax-xmin) / FLOAT(NrPixels);
  133.     dy := (ymax-ymin) / FLOAT(NrPixels);
  134.     y := ymin;
  135.     FOR iy := 0 TO NrPixels-1 DO
  136.       x := xmin;
  137.       FOR ix := 0 TO NrPixels-1 DO
  138.         this := iteration(p,q,x,y);
  139.         SetWindow(w);
  140.         Colour(this);
  141.         Dot(w,ix,iy);
  142.         ResetWindow;
  143.         x := x + dx;
  144.       END;
  145.       y := y + dy;
  146.       BusyRead(ch); IF (ch = ' ') THEN iy := NrPixels END;
  147.     END;
  148.     RestoreFPEnv;
  149.   END Calculate;
  150.  
  151.  
  152.   PROCEDURE ShowFunction(v: Window);
  153.   VAR ix,iy: INTEGER; x,y,y0,fmin,fmax,ry: REAL;
  154.       a: ARRAY [0..NrPixels] OF REAL;
  155.   BEGIN
  156.     SaveFPEnv;
  157.     ClearFPEnv;
  158.     fmin := 0.0; fmax := 0.0;
  159.     FOR ix := 0 TO NrPixels DO
  160.       x := xmin + FLOAT(ix) * ( (xmax-xmin) / FLOAT(NrPixels) );
  161.       y := f(x); a[ix] := y;
  162.       IF y > fmax THEN fmax := y
  163.       ELSIF y < fmin THEN fmin := y
  164.       END;
  165.     END;
  166.     FOR ix := 0 TO NrPixels DO a[ix] := a[ix] - fmin END;
  167.     FOR iy := 0 TO NrPixels DO
  168.       SetWindow(v);
  169.       ForeColour(zyan);
  170.       Dot(v, NrPixels DIV 2, iy);
  171.       ResetWindow;
  172.     END;
  173.     ry := ABS(fmax-fmin);
  174.     y0 := ABS(fmin) / ry * FLOAT(NrPixels);
  175.     iy := TRUNC(y0);
  176.     FOR ix := 0 TO NrPixels DO
  177.       SetWindow(v);
  178.       ForeColour(zyan);
  179.       Dot(v, ix, iy);
  180.       ResetWindow;
  181.     END;
  182.     FOR ix := 0 TO NrPixels DO
  183.       y := (a[ix] / ry) * FLOAT(NrPixels);
  184.       iy := TRUNC(y);
  185.       SetWindow(v);
  186.       ForeColour(red);
  187.       Dot(v, ix, iy);
  188.       ResetWindow;
  189.     END;
  190.     RestoreFPEnv;
  191.   END ShowFunction;
  192.  
  193.  
  194. BEGIN
  195.  
  196.   max := 1.0E-5;
  197.   maxiter := 255;
  198.   xmax := 4.0; xmin := -4.0;
  199.   ymax := 4.0; ymin := -4.0;
  200.   p := 0.0; q := 0.0;
  201.  
  202.   OpenGraphicWindow(v,40,20,NrPixels+4,NrPixels+20,"Funktion",Clear);
  203.   Clear(v);
  204.   ShowFunction(v);
  205.  
  206.   OpenGraphicWindow(w,330,100,NrPixels+2,NrPixels+20,"Newton",Clear);
  207.   Clear(w);
  208.  
  209.   LOOP
  210.     WriteString ('drawing Newton.'); WriteLn;
  211.     Calculate; ch := 0C;
  212.     WriteString ('zoom requested Y/N:');
  213.     Read(ch); Write(ch); WriteLn;
  214.     IF CAP(ch) # "Y" THEN EXIT END;
  215.     WriteString ('define window !'); WriteLn;
  216.     GetPos(ox,oy);
  217.     ux := xmin + FLOAT(ox)*dx;
  218.     uy := ymin + FLOAT(oy)*dy;
  219.     GetPos(nx,ny);
  220.     deltax := ABS(nx-ox);
  221.     uux := xmin + FLOAT(nx)*dx;
  222.     uuy := uy + FLOAT(deltax)*dy;
  223.     xmin := ux; xmax := uux;
  224.     ymin := uy; ymax := uuy;
  225.  
  226.     SetWindow(w);
  227.     ForeColour(black);
  228.     SetPen(w,ox,oy);        MoveTo(w,nx,oy);
  229.     SetPen(w,nx,oy);        MoveTo(w,nx,oy+deltax);
  230.     SetPen(w,nx,oy+deltax); MoveTo(w,ox,oy+deltax);
  231.     SetPen(w,ox,oy+deltax); MoveTo(w,ox,oy);
  232.     ResetWindow;
  233.  
  234.     GetPos(ox,oy); (* wait for mouse click *)
  235.     Clear(w);
  236.   END;
  237.  
  238.   CloseGraphicWindow(w);
  239.   CloseGraphicWindow(v);
  240.  
  241. END Newton.
  242.